clc;clear;clearvars;
% 随机生成10个数据
num_initial = 10;
num_vari = 30;
% 搜索区间
upper_bound = 5.12;
lower_bound = -5.12;
iter = 20000;
varphi = 1.5;
% 随机生成10个数据,并获得其评估值
sample_x = lhsdesign(num_initial, num_vari).*(upper_bound - lower_bound) + lower_bound.*ones(num_initial, num_vari);
sample_y = Rastrigin(sample_x);
Fmin = zeros(iter, 1);
aver_Fmin = zeros(iter, 1);

for n = 1 : 100
    k = 1;
    % 初始化一些参数
    X = sample_x;
    Y = sample_y;
    V = sample_x;
    [fmin, ~] = min(Y);
    fprintf("n: %.4f\n", n);
    UV = V;
    UX = X;
    UY = Y;
    for i = 1 : iter
        PV = [];
        PX = [];
        PY = [];
        while ~isempty(UX)
            m = size(UX, 1);
            k1 = randi(m);
            k2 = randi(m);
            while k2 == k1
                k2 = randi(m);
            end
            v1 = UV(k1, :);
            x1 = UX(k1, :);
            y1 = UY(k1, :);
            v2 = UV(k2, :);
            x2 = UX(k2, :);
            y2 = UY(k2, :);
            r1 = rand(1, num_vari);
            r2 = rand(1, num_vari);
            r3 = rand(1, num_vari);
            % loss就要更新
            if y1 <= y2
                PV = [PV; v1];
                PX = [PX; x1];
                PY = [PY; y1];
                v2 = r1 .* v2 + r2 .* (x1 - x2) + varphi * r3 .* ((x1 + x2) ./ 2 - x2);
                x2 = x2 + v2;
                x2(x2 > upper_bound) = upper_bound;
                x2(x2 < lower_bound) = lower_bound;
                y2 = Rastrigin(x2);
                PV = [PV; v2];
                PX = [PX; x2];
                PY = [PY; y2];
            else
                PV = [PV; v2];
                PX = [PX; x2];
                PY = [PY; y2];
                v1 = r1 .* v1 + r2 .* (x2 - x1) + varphi * r3 .* ((x1 + x2) ./ 2 - x1);
                v1(v1 > upper_bound) = upper_bound;
                v1(v1 < lower_bound) = lower_bound;
                x1 = x1 + v1;
                x1(x1 > upper_bound) = upper_bound;
                x1(x1 < lower_bound) = lower_bound;
                y1 = Rastrigin(x1);
                PV = [PV; v1];
                PX = [PX; x1];
                PY = [PY; y1];
            end
            % 为了删除防止下标发生变化
            if k1 > k2
                UV(k1, :) = [];
                UX(k1, :) = [];
                UY(k1, :) = [];
                UV(k2, :) = [];
                UX(k2, :) = [];
                UY(k2, :) = [];
            else
                UV(k2, :) = [];
                UX(k2, :) = [];
                UY(k2, :) = [];
                UV(k1, :) = [];
                UX(k1, :) = [];
                UY(k1, :) = [];
            end
        end
        UV = PV;
        UX = PX;
        UY = PY;
        % 更新所有个体最佳位置
        [fmin, ~] = min(UY);
        % fprintf("iter %d fmin: %.4f\n", i, fmin);
        Fmin(k, 1) = fmin;
        k = k +1;
    end
    aver_Fmin = aver_Fmin + Fmin;
end
aver_Fmin = aver_Fmin ./ 100;
% disp(pbestx(gbest, :));
plot(aver_Fmin);